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The  maximum  speed  with  which  information  can  propagate  in  a  quan¬ 
tum  many-body  system  directly  affects  how  quickly  disparate  parts  of 
the  system  can  become  correlated1-4  and  how  difficult  the  system  will  be 
to  describe  numerically5.  For  systems  with  only  short-range  interactions, 
Lieb  and  Robinson  derived  a  constant-velocity  bound  that  limits  corre¬ 
lations  to  within  a  linear  effective  Tight  cone’6.  However,  little  is  known 
about  the  propagation  speed  in  systems  with  long-range  interactions, 
because  analytic  solutions  rarely  exist  and  because  the  best  long-range 
bound7  is  too  loose  to  accurately  describe  the  relevant  dynamical  time- 
scales  for  any  known  spin  model.  Here  we  apply  a  variable-range  Ising 
spin  chain  Hamiltonian  and  a  variable-range  XY  spin  chain  Hamiltonian 
to  a  far-from-equilibrium  quantum  many-body  system  and  observe  its 
time  evolution.  For  several  different  interaction  ranges,  we  determine 
the  spatial  and  time-dependent  correlations,  extract  the  shape  of  the  light 
cone  and  measure  the  velocity  with  which  correlations  propagate  through 
the  system.  This  work  opens  the  possibility  for  studying  a  wide  range  of 
many-body  dynamics  in  quantum  systems  that  are  otherwise  intractable. 

Lieb-Robinson  bounds6  have  strongly  influenced  our  understanding  of 
locally  interacting,  quantum  many-body  systems.  They  restrict  the  many- 
body  dynamics  to  a  well-defined  causal  region  outside  of  which  correlations 
are  exponentially  suppressed8,  analogous  to  causal  light  cones  that  arise  in 
relativistic  theories.  These  bounds  have  enabled  a  number  of  important 
proofs  in  condensed-matter  physics5,7,9-11,  and  also  constrain  the  timescales 
on  which  quantum  systems  might  thermalize12-14  and  the  maximum  speed 
that  information  can  be  sent  through  a  quantum  channel15.  Recent  experi¬ 
mental  work  has  observed  linear  (that  is,  Lieb-Robinson-like)  correlation 
growth  over  six  sites  in  a  one-dimensional  quantum  gas16. 

When  interactions  in  a  quantum  system  are  long  range,  the  speed  with 
which  correlations  build  up  between  distant  particles  is  no  longer  guaranteed 
to  obey  the  Lieb-Robinson  prediction.  Indeed,  for  sufficiently  long-range 
interactions,  the  notion  of  locality  is  expected  to  breakdown  completely17. 
Inapplicability  of  the  Lieb-Robinson  bound  means  that  comparatively 
little  can  be  predicted  about  the  growth  and  propagation  of  correlations 
in  long-range-interacting  systems,  although  there  have  been  several 
recent  theoretical  and  numerical  advances2,3,7,17-20. 

Here  we  report  direct  measurements  of  the  shape  of  the  causal  region 
and  the  speed  at  which  correlations  propagate  in  an  Ising  spin  chain  and 
a  newly  implemented  XT  spin  chain.  The  experiment  is  effectively  deco¬ 
herence  free  and  serves  as  an  initial  probe  of  the  many-body  dynamics  of 
isolated  quantum  systems.  Within  this  broad  experimental  framework, 
studies  of  entanglement  growth21,  thermalization12,14  or  other  dynamical 
processes— with  or  without  controlled  decoherence— can  be  realized.  Scal¬ 
ing  such  quantum  simulations  to  larger  system  sizes  is  straightforward 
(Methods),  unlike  ground-state  or  equilibrium  studies  that  typically  must 
consider  diabatic  effects22,23. 

To  induce  the  spread  of  correlations,  we  perform  a  global  quench  by 
suddenly  switching  on  the  spin-spin  couplings  across  the  entire  chain  and 
allowing  the  system  to  evolve  coherently.  The  dynamics  following  a  global 
quench  can  be  highly  non-intuitive;  one  picture  is  that  entangled  quasi¬ 
particles  created  at  each  site  propagate  outwards,  correlating  distant  parts 
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of  the  system  through  multiple  interference  pathways13.  This  process  differs 
substantially  from  local  quenches21,  where  a  single  site  emits  quasiparticles 
that  may  travel  ballistically3,13,  resulting  in  a  different  causal  region  and  prop¬ 
agation  speed  than  in  a  global  quench  (even  for  the  same  spin  model). 

The  effective  spin- 1/2  system  is  encoded  into  the  2Si/2|E—  0,  mF  —  0) 
and  |  F=  1,  mF  =  0)  hyperfine  clock’  states  of  trapped 171  Yb+  ions,  denoted 
|  [)z  and  1 1 )z,  respectively24.  We  initialize  a  chain  of  1 1  ions  by  optically  pump¬ 
ing  to  the  product  state  1 .  .)z  (Fig.  1).  At  time  t  =  0,  we  quench  the 
system  by  applying  phonon-mediated,  laser- induced  forces25-27  to  yield  an 
Ising  or  XY  model  Hamiltonian  (Methods) 

Hising  =  / i,j  rfGj  ( 1 ) 

i<j 

(2) 

i<j 

where  o\  (y  =  x,y,  z)  is  the  Pauli  matrix  acting  on  the  ith  spin,  (in  cyclic 
frequency)  is  the  coupling  strength  between  spins  i  and  j,  and  we  use  units 
in  which  Planck’s  constant  equals  1.  For  both  model  Hamiltonians,  the 


O  O 


Figure  1  |  Sketch  of  experimental  protocol.  Step  (1):  the  experiment  is 
initialized  by  optically  pumping  all  11  spins  to  the  state  |  [)z.  Step  (2):  after 
initialization,  the  system  is  quenched  by  applying  laser-induced  forces  on  the 
ions,  yielding  an  effective  Ising  or  XT  spin  chain  (see  text  for  details).  Step  (3): 
after  allowing  dynamical  evolution  of  the  system,  the  projection  of  each  spin 
along  the  2  direction  is  imaged  onto  a  charge-coupled  device  (CCD)  camera. 
Such  measurements  allow  us  to  construct  any  possible  correlation  function  Cy 
along  z. 
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spin-spin  interaction  matrix  /zy  contains  tunable,  long-range  couplings 
that  fall  off  approximately  algebraically  as  oc  ll\i—j\a  (ref.  26).  We 
vary  the  interaction  range  a  by  adjusting  a  combination  of  trap  and  laser 
parameters22  (Methods),  choosing  a  ~  0.63,  0.83,  1.00  or  1.19  for  these 
experiments. 

After  quenching  to  the  Ising  or  XT  model  with  our  chosen  value  of  a, 
we  allow  coherent  evolution  for  various  lengths  of  time  before  resolving 
the  spin  state  of  each  ion  using  a  charge-coupled  device  camera.  The  exper¬ 
iments  at  each  time  step  are  repeated  4,000  times  to  collect  statistics.  To 
observe  the  build-up  of  correlations,  we  use  the  measured  spin  states  to 
construct  the  connected  correlation  function 

Cy(0  =  (*?('  (3) 

between  any  pair  of  ions  at  any  time.  Because  the  system  is  initially  in  a 
product  state,  Q>;( 0)  =  0  everywhere.  As  the  system  evolves  away  from  a 
product  state,  evaluating  equation  (3)  at  all  points  in  space  and  time  pro¬ 
vides  the  shape  of  the  light-cone  boundary  and  the  correlation  propagation 
velocity  for  our  long-range  spin  models. 

Figure  2  shows  the  results  of  globally  quenching  the  system  to  a  long- 
range  Ising  model  for  four  different  interaction  ranges.  In  each  case,  we 
extract  the  light-cone  boundary  by  measuring  the  time  it  takes  a  cor¬ 
relation  of  fixed  amplitude  (here  Qj  =  0.04—0. ICH3*,  where  C™**  is 
the  largest  connected  correlation  between  two  ions)  to  travel  an  ion-ion 
separation  distance  r.  For  strongly  long-range  interactions  (a  <  1),  we  observe 
accelerating  information  transfer  through  the  chain.  This  fast  propagation 
of  correlations  is  not  surprising,  because  even  the  direct  long-range  coup¬ 
ling  between  distant  spins  produces  correlations  in  a  time  t  oc  1//^  ~  ra. 
However,  increasing  propagation  velocities  quickly  surpass  the  Lieb-Robinson 
velocity  for  a  system  with  equivalent  nearest-neighbour-only  interactions, 
vLR  =  12e/max,  where  e  is  Euler’s  number  and  /max  is  the  maximum  Ising 
coupling  strength  for  a  given  spin-spin  coupling  matrix  (Fig.  2c,  f,  i).  This 
serves  as  experimental  confirmation  that  predictions  based  on  the  Lieb- 
Robinson  result— including  those  that  bound  the  growth  of  entanglement 
or  set  thermalization  timescales—  are  no  longer  applicable  when  interac¬ 
tions  are  sufficiently  long  range. 


For  the  specific  case  of  the  pure  Ising  model,  the  correlations  at  any  time 
can  be  predicted  by  an  exact  analytic  solution18,28: 

C/,;(f)=I  n  COS  [2  (/yt +/;,*)  f] 

+  \  n  cos[2(jiJc-Jj'k)t]  (4) 

Zk^i,j 

—  IT  cos[2//  kt]  II  cos  \2L  kt] 

k^i  ’  sk*j  1  J'  J 

In  equation  (4),  correlations  can  only  build  up  between  sites  i  and  j  that 
are  coupled  either  directly  or  through  a  single  intermediate  spin  k;  pro¬ 
cesses  which  couple  through  more  than  one  intermediate  site  are  pro¬ 
hibited.  For  instance,  if  the  couplings  are  nearest-neighbour-only  then 
Ciyj(t)  =  0  for  all  \i  —  j\  >  2.  This  property  holds  for  any  commuting 
Hamiltonian  (Methods)  and  explains  why  the  spatial  correlations  shown 
in  Fig.  2  become  weaker  for  shorter-range  systems. 

The  products  of  cosines  in  equation  (4)  with  many  different  oscillation 
frequencies  result  in  the  observed  decay  of  correlations  when  t  >0.1/  /max. 
At  later  times,  rephasing  of  these  oscillations  creates  revivals  in  the  spin- 
spin  correlation.  One  such  partial  revival  occurs  at  t=  2.44//max  for  a  = 
0.63  (Extended  Data  Fig.  1),  verifying  that  our  system  remains  coherent  on 
a  timescale  much  longer  than  that  which  determines  the  light-cone  boundary. 

We  repeat  the  quench  experiments  for  an  XT  model  Hamiltonian 
using  the  same  set  of  interaction  ranges  a  (Fig.  3).  Dynamical  evolution 
and  the  spread  of  correlations  in  long- range- interacting  XT  models  are 
much  more  complex  than  in  the  Ising  case  because  the  Hamiltonian  con¬ 
tains  non-commuting  terms.  As  a  result,  there  exists  no  exact  analytic  solu¬ 
tion  comparable  to  equation  (4). 

Compared  with  the  correlations  observed  for  the  Ising  Hamiltonian, 
correlations  in  the  XY  model  are  much  stronger  at  longer  distances  (for 
example,  compare  Fig.  2j  with  Fig.  3j).  Processes  coupling  through  mul¬ 
tiple  intermediate  sites  (which  were  disallowed  in  the  commuting  Ising 
Hamiltonian)  now  have  a  critical  role  in  building  correlations  between 
distant  spins.  These  processes  may  also  explain  our  observation  of  a  steeper 


Figure  2  |  Measured  quench  dynamics  in  a  long-range  Ising  model. 

a-c,  Spatial  and  time- dependent  correlations  (a),  extracted  light-cone 
boundary  (b)  and  correlation  propagation  velocity  (c)  following  a  global 
quench  of  a  long-range  Ising  model  with  a  =  0.63.  The  curvature  of  the 
boundary  shows  an  increasing  propagation  velocity  (b),  quickly  exceeding  the 
short-range  Lieb-Robinson  velocity  bound,  vLR  (c).  Solid  lines  give  a  power-law 
fit  to  the  data,  which  slightly  depends  on  the  choice  of  fixed  contour  Czy. 
d-1,  Complementary  plots  for  a  =  0.83  (d-f),  a  =  1.00  (g-i)  and  a  =  1.19  (j-1). 
As  the  range  of  the  interactions  decreases,  correlations  do  not 


propagate  as  far  or  as  quickly  through  the  chain;  the  short-range  velocity 
bound  vLR  is  not  exceeded  for  our  shortest- range  interaction,  m,  n,  Nearest- 
neighbour  (m)  and  tenth-nearest-neighbour  (n)  correlations  for  our  shortest- 
and  longest- range  interactions  show  excellent  agreement  with  the 
decoherence-free  exact  solution  (with  no  adjustable  parameters) 
from  equation  (4)  (solid).  The  dashed  blue  curves  show  an  improved 
long-range  bound  valid  for  any  commuting  Hamiltonian  (Methods).  Error 
bars,  1  s.d. 


10  JULY  2014  |  VOL  511  I  NATURE  |  199 


©2014  Macmillan  Publishers  Limited.  All  rights  reserved 


RESEARCH 


LETTER 


Ion  separation,  r  Time  (1/Jmax) 


0.2 


0.00 


I 


1  4  7 

Ion  separation,  r 


0  0.08  0.16 
Time(1/JmJ 


j 


a; 

E 

F 


Ion  separation,  r 


Figure  3  |  Measured  quench  dynamics  in  a  long-range  XY  model.  Global 
quench  of  a  long-range  XY  model  with  four  different  interaction  ranges. 
a-1,  Panel  descriptions  match  those  in  Fig.  2.  In  each  case,  when  compared 
with  the  Ising  model,  correlations  between  distant  sites  in  the  XY  model  are 
stronger  and  build  up  more  quickly.  For  the  shortest-range  interaction  (j— 1), 
we  observe  a  faster- than -linear  growth  of  the  light-cone  boundary, 


despite  having  a  >  1;  no  known  analytic  theory  predicts  this  effect, 
m,  n,  Measured  nearest-neighbour  and  tenth-nearest-neighbour  correlations 
closely  match  the  numerical  solution  found  by  evolving  the  Schrodinger 
equation  of  an  XY  model  (equation  (2))  with  no  free  parameters  and  no 
decoherence.  Error  bars,  1  s.d. 


power-law  scaling  of  the  light-cone  boundary  in  the  XY  model.  However, 
without  an  exact  solution  there  is  no  a-priori  reason  to  assume  a  power-law 
light-cone  edge  (used  for  the  fits  in  Fig.  3);  deviations  from  power-law 
behaviour  might  reveal  themselves  for  larger  system  sizes. 

An  important  observation  in  Fig.  3j— 1  is  that  of  faster-than-linear  light- 
cone  growth  for  our  shortest-range  interaction,  with  a  =  1.19.  Although 
faster-than-linear  growth  is  expected  for  a  <  1  (see  discussion  of  Ising 
model),  there  is  no  consensus  on  whether  such  behaviour  is  generically 
expected  for  a  >  1 .  Our  experimental  observation  has  prompted  us  to  nu¬ 
merically  check  the  light-cone  shape  fora  =  1 . 1 9;  we  find  that  faster-than- 
linear  scaling  persists  in  systems  of  up  to  22  spins  before  our  calculations 
break  down  (Extended  Data  Fig.  2). 

Whether  such  scaling  continues  beyond  —30  spins  is  a  question  that, 
at  present,  quantum  simulators  are  best  positioned  to  answer.  In  Figs  2m, 
n  and  3m,  n,  the  excellent  agreement  between  data  and  theory  demon¬ 
strates  that  experiments  produce  the  correct  results  in  a  regime  still  solv¬ 
able  by  classical  computers.  For  larger  systems,  where  numerical  evolution 
of  the  Schrodinger  equation  fails,  the  quality  of  quantum  simulations 
could  still  be  benchmarked  against  the  exact  Ising  solution  of  equation  (4). 
Finding  close  agreement  in  the  Ising  case  would  then  build  confidence  in 
an  XY  model  simulation,  which  cannot  be  validated  by  any  other  known 
method. 

For  the  XY  model,  we  additionally  study  the  spatial  decay  of  correla¬ 
tions  outside  the  light-cone  boundary.  The  data  (Fig.  4)  is  well  described 
by  fits  to  exponentially  decaying  functions.  Recent  theoretical  work20 
predicts  an  initial  decay  of  spatial  correlations  bounded  by  an  expo¬ 
nential,  followed  by  a  power-law  decay;  we  speculate  that  much  larger 
system  sizes  and  several  hundred  thousand  repetitions  of  each  data  point 
(to  reduce  the  shot- noise  uncertainty  sufficiently)  would  be  necessary  to 
see  this  effect. 

A  perturbative  treatment  of  time  evolution  under  the  XT  Hamiltonian 
yields  the  short-time  approximation  for  the  correlation  function  C^(t)  ~ 

( 'Jijt )2.  These  values  are  plotted  as  dashed  lines  along  with  the  data  in 
Fig.  4.  Although  the  perturbative  result  matches  the  data  early  on,  it  fails 
to  describe  the  dynamics  at  longer  evolution  times.  The  discrepancies 
indicate  that  the  light-cone  shapes  observed  in  the  XY  model  are  fun¬ 
damentally  non-perturbative;  rather,  they  result  from  the  build-up  of 
correlations  through  multiple  intermediate  sites  and  cannot  be  described 
by  any  known  analytical  method. 


Figure  4  |  Correlations  and  dynamics  beyond  the  perturbative  regime. 

Decay  of  spatial  correlations  outside  the  light-cone  boundaries  for  a  long-range  XT 
model  with  a  =  0.63, 0.83, 1.00  or  1.19.  The  hatched  region  indicates  the  area  inside 
the  light-cone  boundary  Cy  =  0.15.  The  data  corresponds  to  times  indicated  by 
tickmarks  on  the  left  axis.  Solid  lines  give  an  exponential  fit  to  the  data  and  dashed 
lines  show  the  predictions  from  a  perturbative  calculation.  Perturbation  theory 
does  not  accurately  describe  the  dynamics  at  later  times.  Associated  data  and 
theoretical  results  are  similarly  coloured  to  guide  the  eye.  Error  bars,  1  s.d. 
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We  have  presented  experimental  observations  of  the  causal  region 
and  propagation  velocities  for  correlations  following  global  quenches 
in  Ising  and  XT spin  models.  The  long-range  interactions  in  our  system 
lead  to  a  breakdown  of  the  locality  associated  with  Lieb-Robinson  bounds, 
and  dynamical  evolution  in  the  XY  model  leads  to  results  that  cannot  be 
described  by  analytic  or  perturbative  theory.  Our  work  demonstrates  that 
quantum  simulators  with  only  a  few  tens  of  spins  can  be  an  important  tool 
for  investigating  and  enriching  our  understanding  of  dynamics  in  com¬ 
plex  many-body  systems. 


METHODS  SUMMARY 

We  generate  spin-spin  interactions  by  applying  spin-dependent  optical  dipole 
forces  to  ions  confined  in  a  three-layer  linear  Paul  trap  with  a  4.8  MHz  radial  fre¬ 
quency.  Two  off-resonance  laser  beams  with  a  wavevector  difference  A k  along  a 
principal  axis  of  transverse  motion  globally  address  the  ions  and  drive  stimulated 
Raman  transitions.  The  two  beams  contain  a  pair  of  beat-note  frequencies  symmet¬ 
rically  detuned  from  the  resonant  transition  at  v0  =  12.642819  GHz  by  a  frequency 
p,  comparable  to  the  transverse  motional  mode  frequencies.  In  the  Lamb-Dicke 
regime,  this  results  in  the  Ising- type  Hamiltonian  in  equation  (l)25,26  with 


N 


jUj=a2cDR  Y 

m=  1 


where  Q  is  the  global  Rabi  frequency,  coR  =  hAYl2M  (/z,  Planck’s  constant  divided  by 
2n)  is  the  recoil  frequency,  bi>m  is  the  normal-mode  matrix29  and  com  are  the  transverse 
mode  frequencies.  The  coupling  profile  may  be  approximated  as  a  power-law  decay 
Ji,j  ~  Jo/\i  —j\y'>  where  in  principle  a  can  be  tuned  between  0  and  3  by  varying  the  laser 
detuning  p  or  the  trap  frequencies  0)m  (refs  22,  26). 

We  implement  a  tunable-range  XY  model  by  adding  an  effective  transverse  mag¬ 
netic  field  B  i^i  to  the  pure  Ising  Hamiltonian  with  an  additional  laser  beat-note 
frequency  at  v0.  In  the  limit  B^>J,  processes  governed  by  the  o*  a3-  coupling  which  flip 
two  spins  along  y  (for  example  o+o+,  where  here  o±  =  rf±  icY)  are  energetically 
forbidden,  leaving  only  the  energy-conserving  flip-flop  terms  ( o+o~  +  o~o+).  At 
times  t=  n/B  (with  integer  n),  the  dynamics  of  the  transverse  field  rephases  and  leaves 
only  the  pure  XY  Hamiltonian  of  equation  (2). 

In  the  limit  B  >  r\mQ,  where  rjm  —  A k^/h/2Mcom,  phonon  contributions  from  the 
large,  non-commuting  transverse  field  can  lead  to  unwanted  spin-motion  entangle¬ 
ment  at  the  end  of  an  experiment30.  Therefore,  this  method  of  generating  an  XY 
model  requires  the  hierarchy  J^B<^pmQ  for  all  m.  For  our  typical  trap  parameters, 
/max  ~  400  Hz,  B  ~  4  kHz  and  rjm Q  ~  20  kHz. 


Online  Content  Methods,  along  with  any  additional  Extended  Data  display  items 
and  Source  Data,  are  available  in  the  online  version  of  the  paper;  references  unique 
to  these  sections  appear  only  in  the  online  paper. 
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Generating  spin-spin  couplings.  We  generate  spin-spin  interactions  by  applying 
spin- dependent  optical  dipole  forces  to  ions  confined  in  a  three-layer  linear  Paul  trap 
with  a  4.8  MHz  radial  frequency.  Two  off-resonance  laser  beams  with  a  wavevector 
difference  A  k  along  a  principal  axis  of  transverse  motion  globally  address  the  ions  and 
drive  stimulated  Raman  transitions.  The  two  beams  contain  a  pair  of  beat-note  fre¬ 
quencies  symmetrically  detuned  from  the  resonant  transition  at  v0  =  12.642819  GHz 
by  a  frequency  /i,  comparable  to  the  transverse  motional  mode  frequencies.  In  the 
Lamb-Dicke  regime,  this  results  in  the  Ising-type  Hamiltonian  in  equation  (l)25,26  with 


2  A  b,mbjm 

jij=Q2(ORy2^^ 


(5) 


where  Q  is  the  global  Rabi  frequency,  coR  =  liAk2/2M  is  the  recoil  frequency,  b,„,  is  the 
normal-mode  matrix29  and  (Dm  are  the  transverse  mode  frequencies.  The  coupling 
profile  may  be  approximated  as  a  power-law  decay  JUj  ~  JJ  \  i  —  j\ a,  where  in  principle  a 
can  be  tuned  between  0  and  3  by  varying  the  laser  detuning  /r  or  the  trap  frequencies 
c om  (refs  22,  26). 

We  implement  a  tunable-range  XT  model  by  adding  an  effective  transverse  mag¬ 
netic  field  B  1°  the  pure  Ising  Hamiltonian  with  an  additional  laser  beat-note 
frequency  at  v0.  In  the  limit  processes  in  the  <7? <7  j  coupling  which  flip  two  spins 

along  y  (for  example  o+o+,  where  here  a±  =  o*  ±  ia  )  are  energetically  forbidden, 
leaving  only  the  energy-conserving  flip-flop  terms  {o+o~  +  oo+).  At  times  t  =  n/B 
(with  integer  n),  the  dynamics  of  the  transverse  field  rephases  and  leaves  only  the  pure 
XT  Hamiltonian  of  equation  (2). 

In  the  limit  B  >  r\mQ ,  where  rjm  =  Ak^h/2Mcom,  phonon  contributions  from  the 
large,  non-commuting  transverse  field  can  lead  to  unwanted  spin-motion  entangle¬ 
ment  at  the  end  of  an  experiment30.  Therefore,  this  method  of  generating  an  XT  model 
requires  the  hierarchy  J<^B<^rjm£2  for  all  m.  For  our  typical  trap  parameters,  /max  ~ 
400  Hz,  B  ~  4  kHz  and  rjmQ  «  20  kHz. 

State  detection  and  readout.  After  quenching  to,  and  allowing  time  evolution  under, 
our  spin  Hamiltonian,  we  measure  the  spin  projections  of  each  ion  along  the  z  direc¬ 
tion  of  the  Bloch  sphere.  For  3  ms,  we  expose  the  ions  to  a  laser  beam  that  addresses 
the  cycling  transition  2Si/2  |  F  =  1)  to  2P1/2 1 F  =  0).  Ions  fluoresce  only  if  they  are  in  the 
state  1 1 )z-  This  fluorescence  is  collected  through  an  objective  with  a  numerical  aper¬ 
ture  of  0.23  and  imaged  using  an  intensified  CCD  camera  with  single-site  resolution. 

To  discriminate  between  ‘bright’  and  ‘dark’  states  ( |  ])z  and  |  [)z,  respectively),  we 
begin  by  calibrating  the  camera  with  1,000  cycles  each  of  all-bright  and  all-dark  states. 
For  the  bright  states,  the  projection  of  the  two-dimensional  CCD  image  onto  a  one¬ 
dimensional  (ID)  row  gives  a  profile  comprised  of  Gaussians  at  each  ion  location.  We 
perform  fits  to  locate  the  centre  and  fluorescence  width  of  each  ion  on  our  CCD. 

We  achieve  single- shot  discrimination  of  individual  ion  states  in  the  experimental 
data  by  fitting  the  captured  ID  profile  to  a  series  of  Gaussians  with  calibrated  widths 
and  positions  but  freely  varying  amplitudes.  The  extracted  amplitudes  for  each  ion  are 
then  compared  with  a  threshold  found  by  Monte  Carlo  simulation  to  determine  whether 
the  measured  state  was  bright  or  dark.  Our  discrimination  protocol  also  gives  an  estimate 
of  the  detection  error  (for  example,  misidentifying  a  bright  ion  as  dark),  which  is  typically 
of  order  5%.  Corrected  state  probabilities  (along  with  their  respective  errors)  are  found 
following  the  method  outlined  in  ref.  31,  which  also  takes  into  account  errors  due  to 
quantum  projection  noise. 

Quantum  coherence  and  scaling.  For  the  data  presented  in  Figs  2  and  3  for  the  Ising 
and  XT  models,  correlations  propagate  across  the  entire  chain  in  a  time  t «  0.1  //max. 
During  this  time,  we  find  excellent  agreement  between  the  data  and  decoherence-free 
numerical  simulations.  This  indicates  that  the  experiment  remains  coherent  on  the  time- 
scales  needed  to  measure  the  light-cone  shape  and  correlation  propagation  velocity. 

To  investigate  coherences  at  longer  times,  we  look  for  partial  revivals  of  correla¬ 
tions  as  predicted  by  equation  (4),  with  a  =  0.63.  (Other  values  of  a  are  not  predicted 
to  show  partial  revivals  of  measurable  size  within  experimentally  accessible  timescales.) 
Extended  Data  Fig.  1  shows  evidence  that  the  system  remains  coherent  until  at  least 
t  ~  2.5 //max,  substantially  longer  than  is  needed  to  extract  the  light-cone  boundary  for 
an  11 -spin  system. 

A  natural  question  to  ask  is  how  long  it  takes  to  observe  the  full  light-cone  boundary 
in  an  IV-spin  system.  We  may  then  estimate  the  potential  for  scaling  to  larger  numbers 
while  keeping  the  demonstrated  coherence  time  and  all  other  experimental  para¬ 
meters  fixed.  In  the  worst  case,  an  Ising  model  with  a  =  1.19,  the  light-cone  boundary 
spreads  across  five  sites  in  a  time  t  =  0.0 7//max,  and  grows  as  r  oc  t0'97.  If  correlations 
were  to  continue  spreading  at  this  rate,  they  would  reach  nearly  100  sites  away  within 
our  demonstrated  t  =  2.5//max  coherence  time.  If  one  is  instead  interested  in  seeing 
whether  faster-than-linear  growth  persists  at  the  same  rate  for  an  XT  model  with  a  = 
1.19,  correlations  could  potentially  spread  over  700  sites  within  t  =  2.5/ /max.  Although 
technical  challenges  certainly  must  be  addressed  before  scaling  to  these  large  sizes,  we 
note  that  the  state  initialization,  evolution,  and  measurement  procedures  implemen¬ 
ted  here  would  remain  unchanged  even  for  hundreds  of  ions. 


Lieb-Robinson  velocity  for  nearest-neighbour  interactions.  Here  we  justify  our 
claim  that  the  Lieb-Robinson  velocity6  for  the  spread  of  correlation  functions  from 
an  initial  product  state,  evolving  under  a  ID  spin  Hamiltonian  with  only  nearest- 
neighbour  interactions,  is  bounded  above  by  vLR  =  12c/.  In  particular,  we  consider 
a  Hamiltonian 

»=E**/+» 

j 

with  interaction  strength  |  \hjtj+i\  |  =  /.  Note  that  both  the  Ising  and  XT  Hamiltonians 
defined  in  the  manuscript  satisfy  these  assumptions  in  the  a— >  oo  limit,  where  Jy  = 

jSjj+i,  as  can  easily  be  checked  by  calculating  1 1  °]  1 1  —  1 1  ^  1 1  j ^  —  1-  F°r 

operators  evolving  in  the  Heisenberg  picture  under  H  (such  that  A(t)  =  elHtA(0)e~lHt), 
we  would  like  to  compute  the  connected  correlation  function 

CiJ(t)  =  (Ai(t)Bj(t))c 

where  A{  and  Bj  are  supported  on  sites  i  and/,  respectively. 

A  bound  on  these  correlation  functions  follows  immediately  from  results  in  ref. 
8,  which  relate  a  Lieb-Robinson  bound  on  unequal-time  commutators  to  a  bound 
on  connected  correlation  functions.  In  particular,  for  a  Lieb-Robinson  commutator 
bound  of  the  form 

||[Ai(0,Bi(0)]||sc||A,||||B>||e(rf-^ 

we  have 


Cy(f)<4c||Ai||||B;||e<v,-r/2>/{  (6) 

where  r  is  the  distance  between  the  two  sites  i  and/. 

The  Lieb-Robinson  commutator  bound  for  a  nearest- neighbour  Hamiltonian 
on  a  D-dimensional  square  lattice  is  given  by8 


||[A((t)^(0)]|lS2||Al||||B, 


X 


(2fr(4D-l))fc 

k\ 


which  in  ID  gives 

||  [Aj(f),B,-(0)J|<2||Ai||||Bj||e_r  53^“Tj“^“ 

k  =  r 

<2||A||(||B/||e6e/i“r 

and,  hence,  v  =  6eJ.  The  velocity  bound  for  the  spreading  of  correlations  is  obtained 
by  setting  the  bound  on  Ckft)  (the  right-hand  side  of  equation  (6))  to  a  constant  value 
and  solving  r  =  vLRf,  which  yields  vLR  =  2v  =  12c/. 

Bound  for  commuting  Hamiltonians.  Motivated  by  applications  to  the  Ising  model 
studied  in  the  manuscript,  here  we  derive  a  bound  applicable  to  ID  Hamiltonians 

H=J> 

k<l 

where  [hkb  h^y]  =  0  for  any  k,  /,  kr,  l'.  As  above,  we  are  interested  in  bounding  the 
connected  correlation  function  Cuft),  and  without  loss  of  generality  we  take  i  </.  Lor 
convenience  in  what  follows,  we  define  hkk  =  0,  and  take  hkj  =  hjk  (even  though  only 
one  of  the  two  appears  in  the  Hamiltonian). 

To  compute  A t{t),  let  us  first  define  Ht  =  J2k  hk  as  the  part  of  H  that  (in  general) 
does  not  commute  with  Ab  so  that  A;  (t)  =  elH,tAie~lH,t.  We  can  further  separate  H{ 
into  two  parts  by  choosing  a  site  index  k0  satisfying  i  <  k0  </  and  writing 

H'=  ^  ft* 

k<k0 

H';=  53% 

k>k0 


As  a  result 


iH ' t  JH" t  A  iH' t 


Ai(t)=eiH>teiHitAie 


=  etHS  Af  + 


mA'(t)+fi(t) 

where  A'(i)  =  e^'A^-^'  and 

ll/i(t)ll^2f||A,|[||Hf|[ 
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Similarly,  we  can  define 

k>k0 

k<k0 

and  Aj(t)  =  etHjtAje~lHjti  such  that  Aj(t)=A'j(t)  +fi(t)  and 

l^'W  ||  <2t|h||  l^'l 

In  terms  of  these  newly  defined  quantities,  we  can  write 

Qj(t) = (A'wy(f))c+(/^ (t)Aj(t))c+(A'l(t)fj(t))c 
where  we  note  that  the  second  term  contains  Aj(t)  (rather  than  A'-(t)).  By  inspec¬ 
tion,  (f)Aj  ( t )  ^  =  0.  Using  the  bounds  on  |  \  |  and  |  \  | ,  together  with  the 

inequality  |(AB)C|  <  2| \A\ \ \ \B\\,  we  find  that 

|Cy(f)|<4f||Ai||||A;.||(||H;'||  +  ||^"||) 

Noting  that  \Jki\  =  | \hki\\,  we  then  have 

lQ,W|^iNiwfEw+EW) 

\k>k0  k<k0  ) 

One  can  optimize  the  value  of  k0  to  give  the  tightest  bound.  For  power-law  couplings 
Jki~  \k~l\ -a  (a  >  0)  in  ID,  choosing  k0  right  in  the  middle  of  i  and  j  will  generally 
give  the  tightest  bound. 

Multi-hop  processes  are  forbidden  for  commuting  Hamiltonians.  Here  we  prove 
the  claim  that,  given  an  initial  product  state  evolving  under  a  commuting  Hamiltonian, 
distant  spins  can  only  become  correlated  if  they  are  either  directly  coupled  or  if  they 
share  an  intermediate  spin  to  which  they  both  couple;  multi-hop  processes  (for  example 
site  A  coupling  to  site  D  through  sites  B  and  C)  do  not  occur. 

We  consider  the  time  evolution  of  the  operators  At  and  Ap  residing  on  sites  i  and 
j  of  the  lattice.  As  discussed  in  the  previous  section,  the  time  evolution  of  A{  and  Aj 
can  be  written  as 

Ai{t)  =  eiHitAie~iHit 
Aj  (t)  =  eiHjtAj  e  “ iHjt 

where 

P 

Hj  —  ^2  hjq 


We  can  expand  the  time-evolution  operator  to  obtain 

Ai(t)  =  At  +  ih[HhAi]  ~  f  [Hi, [Hi  A-]]  +  •  •  • 

,2  (7) 

= A,- + it  ^2  [hi  ’M -^£  [V’  [V  >A<]  ]  +  ■  ■  • 

pi  Php2 

It  follows  from  equation  (7)  that  A,-(f)  is  supported  on  (that  is,  can  be  written  in  terms 
of  operators  belonging  to)  site  i  and  any  sitep  for  which  |  \hip\  ]  #  0;  we  denote  the  set 
of  such  points  by  Ab  and  define  an  equivalent  set  Aj  containing  all  sites  supporting 
the  operator  Afjt).  If  | \hy\  |  =  0  and  there  are  no  sites  p  that  simultaneously  satisfy 
IM*°  and  \  \hjp\\  A  0,  then  Aj(AAj  =  0.  In  this  case,  it  is  clear  that  an  initial 
product  state  must  satisfy  ( Ai(t)Aj(t ))  =  (A2(t))(A;(t)),  and  therefore  any  connected 
correlation  function  Cuft)  must  vanish. 

Numerical  solutions.  Because  no  analytic  solution  exists  for  the  XY  model,  exact 
long-time  dynamics  (where  the  perturbative  results  derived  above  break  down)  must 
be  obtained  by  numerical  solution  of  the  Schrodinger  equation.  The  curves  presented 
in  Fig.  3m,  n  are  calculated  using  a  standard  numerical  integration  routine.  With  our 
experimental  spin-spin  couplings  Jj  as  inputs  (equation  (5)),  we  construct  the  full  XY 
Hamiltonian  (equation  (2))  using  sparse  matrices.  After  evolving  the  initial  product 
state  |  t/dO))  under  the  Hamiltonian  HXy for  a  time  t,  we  construct  the  desired  correlation 
functions  by  calculating 


cy(0  =  {^(0|°fcrj|'A(0) 

To  numerically  check  the  light-cone  shape  when  a  =  1.19  in  a  system  of  22  spins, 
we  follow  a  similar  procedure  to  calculate  the  time-evolved  state  |  i/df)).  The  results 
of  this  calculation  are  shown  in  Extended  Data  Fig.  2.  Note  that  faster-than-linear 
growth  of  the  light-cone  boundary  persists  in  this  larger  system  of  22  spins. 
Short-time  perturbation  theory  for  the  XY  model.  Unlike  in  the  Ising  model,  no 
exact  analytic  solution  exists  for  the  XT  model  (even  in  1 D,  owing  to  the  long-range 
couplings).  However,  we  can  nevertheless  expand  the  time-evolution  operator  to 
low  order  and  thereby  recover  the  dynamics  at  short  times.  At  sufficiently  long  times, 
this  perturbative  expansion  (carried  out  here  to  second  order)  becomes  a  poor  approxi¬ 
mation.  This  failure,  which  is  observed  in  the  experimental  dynamics  (Fig.  4),  suggests 
that  the  growth  of  correlations  at  long  distances  is  not  the  result  of  direct  spin-spin 
interactions;  instead  those  correlations  originate  from  the  repeated  propagation  of 
information  through  intermediate  spins. 

We  are  interested  in  the  time  evolution  of  a  connected  correlation  function 
Cij(t)  =  (Ai(t)Aj(t))c  of  observables  At  and  Aj  located  at  different  sites  i  and;.  To 
second  order  in  time,  we  have 

Mt)=A,  +  it[H,Ai]  -  f  [H,[H,Aj]]  +  0(f) 

which  yields 

=  (A‘Ai)c  +  U  {(A‘  lH-Aj]  )c  +  ( [H,A;]Aj)c) 

f2  /  \  (8) 

-  2  (<M«.M]>cU«Aj]]A;->c) 

■—  f2(  [H,Aj]  [H,Aj]  )c  +  0(f3) 

Note  that  in  equation  (8)  we  use  the  notation 

(Ai[HyAj\)c  =  (Aj[HtAj])  -  ( At)([H,Aj ]) 

In  the  experiment,  where  Af  corresponds  to  the  Pauli  spin  operator  0,  the  initial 
state  is  (1)  a  product  state  |  J,  •  •  •  [)z  and  (2)  a  simultaneous  eigenstate  of  each  Ab  As 
a  result  of  (1),  the  connected  correlation  at  t  =  0  vanishes  ((A2A;)C  =  0).  As  a  result 
of  (2),  the  third  and  fourth  lines  in  equation  (8)  vanish.  Therefore,  we  have 

(<7f(0<7;Z(f))c=-«2([H,<7j]  + 

For  the  XT  Hamiltonians  we  find 

[H,ozj]  =  -iJ2Jikhh 

k^i 

and  so 

=t2  E  Mf((^^)A^)(^<7?))  w 

k  ^  i,l  7*;' 

+o{h 

Because  the  initial  state  is  polarized  along  z,  the  only  term  that  has  a  non- zero 
expectation  value  on  the  right-hand  side  of  equation  (9)  is  the  one  with  k  =  j  and 
/  =  i.  Therefore 

=t2j?j(hh)+°(A 
=  (¥)2+o(A 

which  is  the  short-time  result  used  in  the  main  text. 

31.  Shen,C.&  Duan,  L.-M.  Correcting  detection  errors  in  quantum  state  engineering 
through  data  processing.  New  J.  Phys.  14,  053053  (2012). 
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Extended  Data  Figure  1  |  A  long-time  partial  revival  in  the  long-range  Ising  sites  1  and  2  is  evident,  showing  quantum  coherence  at  long  times.  The  black 

model,  a,  Spatial  correlations  measured  at  long  times  after  a  global  quench  of  line  shows  the  exact  solution  predicted  from  equation  (4).  Error  bars,  1  s.d. 

an  Ising  model  with  a  =  0.63.  b,  A  small  partial  revival  in  correlation  between 
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Extended  Data  Figure  2  |  Numeric  calculation  of  XY  model  correlations. 

Calculated  spatial  and  time-dependent  correlations  for  an  N  =  22-spin  XY 
model  with  spin-spin  couplings  Jy  ~  J0/ 1  i  —  j\ 119 ,  found  by  numerically 
evolving  the  Schrodinger  equation. 
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